Testing Magnetospheric Accretion as an Hα Emission Mechanism of Embedded Giant Planets: The Case Study for the Disk Exhibiting Meridional Flow Around HD 163296

When planets accrete the surrounding gas with an accretion rate of ${\dot{M}}_{{\rm{p}}}$, the total accretion luminosity (L_acc) is given as

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{acc}} & = & \displaystyle \frac{{{GM}}_{{\rm{p}}}{\dot{M}}_{{\rm{p}}}}{{R}_{{\rm{p}}}}\\ & \simeq & 1.2\times {10}^{-4}{L}_{\odot }\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)\left(\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{{10}^{-7}{M}_{{\rm{J}}}\,{{\rm{yr}}}^{-1}}\right){\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 1 }$$

where L_⊙ is the solar luminosity, M_p and R_p are the planet mass and radius, and M_J and R_J are Jupiter's mass and radius, respectively. Hereafter, we adopt R_p = 1.5R_J (Table 1); the radius evolution of planets becomes minimal after the initial contraction of planetary envelopes with the Kelvin–Helmholtz timescale ends (e.g., Bodenheimer et al. 2000; Mordasini et al. 2012b). In the subsequent stage, planets undergo the so-called disk-limited gas accretion (e.g., Hasegawa et al. 2019), and rotationally supported circumplanetary disks should emerge due to the conservation of angular momentum. Magnetospheric accretion may come into play at the disk-limited gas accretion stage. In Equation (1), we adopt characteristic values for M_p and ${\dot{M}}_{{\rm{p}}}$ suggested for PDS 70 b/c as an example (Hasegawa et al. 2021 and references herein); these values vary in the following sections.

The presence of circumplanetary disks divides the accretion luminosity into two components: The luminosity from the disks (L_disk = f_in L_acc), and the luminosity from energy that is liberated as the gas of the disk falls onto planets (L_gas = (1 − f_in)L_acc), where f_in is the partition coefficient of the accretion energy that is controlled by the location of the inner disk edge (R_in). When disks are heated predominantly by viscosity, f_in is written as (e.g., Pringle 1981)

$$\begin{eqnarray}&&{f}_{\mathrm{in}}=\displaystyle \frac{3}{2}\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}\right)\left[1-\displaystyle \frac{2}{3}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}\right)}^{1/2}\right].\end{eqnarray} \tag{ 2 }$$

For simplicity, we adopt the above expression for f_in; if R_in = R_p, then f_in = 1/2. In this paper, L_disk and L_gas are referred to as the disk and infall gas luminosities, respectively.

The infall gas luminosity can be further decomposed into two components when planets undergo magnetospheric accretion, that is, R_in > R_p; for this case, infall gas is channeled by the magnetic field lines of accreting planets, and it reaches the planetary surfaces at nearly freefall velocity (see Section 2.5). The gas radiates some energy both when it is in the magnetospheric flow and when it produces a shock at the planetary surfaces. The corresponding luminosity can be written as L_rad = (1 − f_L)L_gas, where f_L is the partition coefficient of L_gas. It is important that not all of the energy of the infall gas radiates away, and hence some of the energy should be thermalized with the atmospheric gas of planets (Marleau et al. 2019). Accordingly, accreting giant planets are heated up by the gas from circumplanetary disks to some extent. We label this luminosity L_therm = f_L L_gas. Aoyama et al. (2020) estimate the value of f_L due to accretion shock at the planetary surfaces and find that f_L ≃ 0.3–0.8.

In summary, the accretion luminosity can be written as

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{acc}} & = & {L}_{\mathrm{disk}}+{L}_{\mathrm{gas}}\\ & = & {L}_{\mathrm{disk}}+{L}_{\mathrm{rad}}+{L}_{\mathrm{therm}},\end{array}\end{eqnarray} \tag{ 3 }$$

where L_disk = f_in L_acc, L_rad = (1 − f_L)(1 − f_in)L_acc, and L _therm = f_L(1 − f_in)L_acc. Figure 1 shows the change of each component of the luminosities as a function of f_in (or R_in) and f_L. As expected, L_gas exceeds L_disk when R_in > R_p. Moreover, L_therm can contribute about up to 70% of L_acc when f_L = 0.8. This indicates that the effective temperature of accreting planets is considerably affected by accretion.

Zoom In Zoom Out Reset image size

In the following calculations, we adopt f_L = 0.5 as it is an intermediate value.

The effective temperature is one important quantity for characterizing the properties of accreting planets. In this work, it becomes the key parameter for estimating planetary magnetic fields. We here compute the effective temperature of accreting giants using L_therm, as discussed above.

The quantification of the effective temperature of young giant planets receives considerable attention in the literature because it may be used as a diagnostic to differentiate their formation mechanisms (e.g., Marley et al. 2007; Spiegel & Burrows 2012). The so-called hot and cold starts (i.e., large-sized planets with high temperatures and small-sized planets with low temperatures) may be realized as the results of two competing planet formation scenarios: gravitational instability and core accretion. The advent of the direct-imaging technique in the search for young giant planets enables measurements of their luminosity, and hence their size and temperature can be estimated (e.g., Marois et al. 2008). This capability thus makes it possible to specify the formation mechanisms of these planets. However, no conclusive remark has yet been made in the literature.

In this work, the effective temperature of accreting giants is computed as

$$\begin{eqnarray}&&{T}_{{\rm{p}},{\rm{e}}}^{4}={T}_{\mathrm{int}}^{4}+\displaystyle \frac{{L}_{\mathrm{therm}}}{4\pi {\sigma }_{\mathrm{SB}}{R}_{{\rm{p}}}^{2}},\end{eqnarray} \tag{ 4 }$$

where T_int is the intrinsic temperature of planets, and σ_SB is the Stefan-Boltzmann constant. The above equation assumes that accreted gas with a luminosity of L_therm is thermalized over the entire surface of the planets. This is most conservative because T_p,e takes the lowest value. An inclusion of L_therm in Equation (4) corresponds to the so-called warm start as some of accreted gas heats the planets. Reliable calculations of T_int require tracking the planet formation histories from the beginning, as was done by Mordasini et al (e.g., 2012a), which is beyond the scope of this work. We therefore assume that T_int = 700 K, following Spiegel & Burrows (2012).

Heating by L_therm leads to the following planet surface temperature:

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{therm}} & \equiv & {\left(\displaystyle \frac{{L}_{\mathrm{therm}}}{4\pi {\sigma }_{\mathrm{SB}}{R}_{{\rm{p}}}^{2}}\right)}^{1/4}\\ & \simeq & 1.5\times {10}^{3}{\rm{K}}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{1/4}{\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{1/2}\right)}^{1/4}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/4}{\left(\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{{10}^{-7}{M}_{{\rm{J}}}\,{{\rm{yr}}}^{-1}}\right)}^{1/4}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-3/4}.\end{array}\end{eqnarray} \tag{ 5 }$$

We use both T_int and T_therm below to explore the strength of the planetary magnetic fields that is produced by these temperatures.

Planetary magnetic fields are generated by dynamo activities operating in electrically conducting interiors, where convective motion occurs. A precise determination of the field strength is hard. To make the problem tractable, we here use a scaling law that is available in the literature and estimate the strength of the planetary magnetic fields.

In principle, the ultimate source of energy to invoke dynamo activities is the thermodynamic energy available in the interior of planets; the energy is converted into magnetic energy, and thermal flux is maintained against ohmic dissipation. Christensen et al. (2009) adopt this principle and derive a scaling law, which is written as

$$\begin{eqnarray}&&\displaystyle \frac{\langle B{\rangle }^{2}}{8\pi }={{cf}}_{\mathrm{ohm}}\langle \rho {\rangle }^{1/3}{\left({Fq}\right)}^{2/3},\end{eqnarray} \tag{ 6 }$$

where 〈B〉 is the mean magnetic field on the dynamo surface, c is the constant of proportionality, f_ohm ≃ 1 is the ratio of ohmic dissipation to total dissipation, 〈ρ〉 is the mean bulk density of planets where the field is generated, and F = 0.35 is the efficiency factor of converting thermal energy into magnetic energy, $q={\sigma }_{\mathrm{SB}}{T}_{{\rm{p}},{\rm{e}}}^{4}$. A value of c ≃ 1.1 is obtained by adopting the typical values of Jupiter (B_{p, s} = 10 G and q = 5.4 × 10³ erg s⁻¹ cm⁻²) and assuming that 〈B〉/B_{p, s} ≃ 7, where B_{p, s} is the magnetic field strength at the planetary surfaces. Remarkably, Christensen et al. (2009) show that the law successfully reproduces the magnetic fields of objects reasonably well, from solar system planets (e.g., Earth and Jupiter) up to rapidly rotating stars such as CTTSs.

We use the above scaling law to compute the magnetic field strength of accreting giants. Combing Equations (4) and (6), we write the magnetic fields of accreting giant planets as

$$\begin{eqnarray}\begin{array}{rcl}{B}_{{\rm{p}},\ {\rm{s}}} & \simeq & \displaystyle \frac{\langle B\rangle }{7}\\ & = & \displaystyle \frac{{6}^{1/6}}{7}{\left(8\pi \right)}^{1/3}{\left({{cf}}_{\mathrm{ohm}}\right)}^{1/2}{\left(F{\sigma }_{\mathrm{SB}}\right)}^{1/3}{M}_{{\rm{p}}}^{1/6}{R}_{{\rm{p}}}^{-1/2}{T}_{{\rm{p}},{\rm{e}}}^{4/3}.\end{array}\end{eqnarray} \tag{ 7 }$$

When two limits are considered for T_p,e, B_{p, s} is rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{B}_{{\rm{p}},\ {\rm{s}}}^{\mathrm{int}} & \simeq & 1.6\times {10}^{2}{\rm{G}}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/6}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-1/2}{\left(\displaystyle \frac{{T}_{\mathrm{int}}}{700{\rm{K}}}\right)}^{4/3}\end{array}\end{eqnarray} \tag{ 8 }$$

when T_p,e ≃ T_int, and

$$\begin{eqnarray}\begin{array}{rcl}{B}_{{\rm{p}},\ {\rm{s}}}^{\mathrm{therm}} & \simeq & 3.3\times {10}^{2}{\rm{G}}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{1/3}{\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{1/2}\right)}^{1/3}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/2}{\left(\displaystyle \frac{{\dot{M}}_{{\rm{p}}}}{{10}^{-7}{M}_{{\rm{J}}}\,{{\rm{yr}}}^{-1}}\right)}^{1/3}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-3/2}\end{array}\end{eqnarray} \tag{ 9 }$$

when T_p,e ≃ T_therm.

These calculations indicate that when T_p,e ≃ T_int, T_int becomes the fundamental parameter for determining B_p,s (Equation (8)); equivalently, earlier formation histories dictate whether magnetospheric accretion occurs. On the other hand, when T_p,e ≃ T_therm, all the key quantities (e.g., T_p,e, B_p,s, and ${\dot{M}}_{{\rm{p}}}$) can be computed self-consistently (see Equations (5) and (9), and also see Equation (13) as discussed below). This essentially means that disk-limited gas accretion can become energetic enough for physical parameters to self-regulate by the corresponding heating; if magnetospheric accretion operates in this gas accretion stage, the resulting observables (e.g., hydrogen emission lines) serve as a direct probe of the stage.

In the following sections, we consider two limiting cases: T_p,e ≃ T_int and T_p,e ≃ T_therm, and we explore under which conditions giant planets undergo magnetospheric accretion and when (observable) hydrogen lines can be emitted.

Magnetospheric accretion is currently a leading hypothesis to explain the observed Hα emission from PDS 70 b/c (e.g., Aoyama & Ikoma 2019; Thanathibodee et al. 2019; Hasegawa et al. 2021). This accretion mode takes action when the magnetic fields of the planets are strong enough to truncate accreting circumplanetary disks (e.g., Zhu 2015; Batygin 2018; Hasegawa et al. 2021). In this section, we determine when this condition is met.

Circumplanetary disks are truncated when the magnetic pressure (${B}_{{\rm{p}}}^{2}/8\pi$) of the host planets exceeds the ram pressure of the accreting disks (Ghosh & Lamb 1979). Mathematically, it is written as

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{{\rm{p}}}^{2}}{8\pi }={f}_{\mathrm{ram}}{\rho }_{\mathrm{ram}}{v}_{\mathrm{Kep}}^{2},\end{eqnarray} \tag{ 10 }$$

where ${v}_{\mathrm{Kep}}=\sqrt{{{GM}}_{{\rm{p}}}/R}$ is the Keplerian velocity around the planets, R is the distance at the disk midplane measured from the planet center, and ${\rho }_{\mathrm{ram}}\sim {\dot{M}}_{{\rm{p}}}/(4\pi {R}^{2}{v}_{\mathrm{Kep}})$ is the ram pressure of the disks. A value of ${f}_{\mathrm{ram}}=1/\sqrt{2}$ is adopted, following Ghosh & Lamb (1979). We also assume that the magnetic field (B_p) of the planets may be represented well as dipole, that is,

$$\begin{eqnarray}&&{B}_{{\rm{p}}}(r)={B}_{{\rm{p}},{\rm{s}}}{\left(R/{R}_{{\rm{p}}}\right)}^{-3}.\end{eqnarray} \tag{ 11 }$$

From Equations (7), (10), and (11), one can derive a relationship between M_p and ${\dot{M}}_{{\rm{p}}}$ for a given value of R. Considering the two limits for T_p,e, ${\dot{M}}_{{\rm{p}}}$ is given as

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{{\rm{p}}}^{\mathrm{int}} & \simeq & 2.4\times {10}^{-8}{M}_{{\rm{J}}}\,{{\rm{yr}}}^{-1}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{-1/6}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{3/2}{\left(\displaystyle \frac{{T}_{\mathrm{int}}}{700{\rm{K}}}\right)}^{8/3}{\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-7/2}\end{array}\end{eqnarray} \tag{ 12 }$$

when T_p,e ≃ T_int, and

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{{\rm{p}}}^{\mathrm{therm}} & \simeq & 2.5\times {10}^{-7}{M}_{{\rm{J}}}\,{{\rm{yr}}}^{-1}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-3}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)}^{2}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{3/2}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-3/2}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-21/2}\end{array}\end{eqnarray} \tag{ 13 }$$

when T_p,e ≃ T_therm. Note that in the above calculations, we set that R = R_in = 4R_p, that is, f_in = 1/4; equivalently, the truncation radius of the disks due to planetary magnetospheres corresponds to their inner edge. Furthermore, ${\dot{M}}_{{\rm{p}}}^{\mathrm{therm}}$ is computed self-consistently, and hence it becomes a function of M_p and R_p (Equation (13)).

We are now in a position to determine under which condition planetary magnetic fields become strong enough to truncate circumplanetary disks. To proceed, we examine all the quantities (T_p,e, B_p,s, ${\dot{M}}_{{\rm{p}}}$, and L_acc) considered so far. As discussed above, T_int is the fundamental parameter when T_p,e ≃ T_int, while these quantities are all computed self-consistently when T_p,e ≃ T_therm.

We first summarize the relevant equations. When T_p,e ≃ T_int, T_int = 700 K, B_p,s is described by Equation (8), ${\dot{M}}_{{\rm{p}}}$ is given by Equation (12), and L_acc is written as

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{acc}}^{\mathrm{int}} & \simeq & 4.4\times {10}^{-5}{L}_{\odot }{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{5/6}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{1/2}{\left(\displaystyle \frac{{T}_{\mathrm{int}}}{700{\rm{K}}}\right)}^{8/3}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-7/2}.\end{array}\end{eqnarray} \tag{ 14 }$$

When T_p,e ≃ T_therm,

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{therm}} & \simeq & 1.7\times {10}^{3}{\rm{K}}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-3/4}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{3/4}\\ & & \times {\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)}^{3/4}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{5/8}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-9/8}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-21/8},\end{array}\end{eqnarray} \tag{ 15 }$$

where Equations (5) and (13) are used,

$$\begin{eqnarray}\begin{array}{rcl}{B}_{{\rm{p}},\ {\rm{s}}}^{\mathrm{therm}} & \simeq & 5.1\times {10}^{2}{\rm{G}}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-1}\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)\\ & & \times \left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right){\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-2}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-21/6},\end{array}\end{eqnarray} \tag{ 16 }$$

where Equations (9) and (13) are used, ${\dot{M}}_{{\rm{p}}}^{\mathrm{therm}}$ is given by Equation (13), and

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{acc}}^{\mathrm{therm}} & \simeq & 4.7\times {10}^{-4}{L}_{\odot }{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-3}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)}^{2}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{5/2}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-5/2}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-21/2},\end{array}\end{eqnarray} \tag{ 17 }$$

where Equations (1) and (13) are used.

We then explore how these quantities behave as a function of M_p for the given values of R = R_in. Figure 2 shows the results. In the plot, two values of R_in are considered: R_in = 2R_p(=3R_J), and 4R_p(=6R_J). It is obvious that when T_p,e ≃ T_int, both T_int and ${B}_{{\rm{p}},{\rm{s}}}^{\mathrm{int}}$ are independent of R_in (see the dashed lines in the two top panels); for ${B}_{{\rm{p}},{\rm{s}}}^{\mathrm{int}}$, it becomes a weak function of M_p. This is simply because the planetary magnetic fields are mainly regulated by the effective temperature (see Equation (7)). When T_p,e ≃ T_therm, both T_therm and ${B}_{{\rm{p}},{\rm{s}}}^{\mathrm{therm}}$ become a decreasing function of R_in and an increasing function of M_p. This is the direct outcome that these solutions are obtained self-consistently; a low value of R_in means a small truncation radius of the disks. This situation tends to occur when the accretion rate is high. To maintain disk truncation against the resulting high ram pressure, the magnetic fields need to be strong, which in turn requires high effective temperatures. For the M_p dependence, it can be understood as follows: massive planets have a deep gravitational potential, which leads to a high ram pressure. In order to prevent the inner disk edge from clashing onto the planetary surfaces, high magnetic fields and hence high effective temperatures are required.

Zoom In Zoom Out Reset image size

The behaviors of ${\dot{M}}_{{\rm{p}}}$ and L_acc are explained in a similar way (the two bottom panels). Note that when T_p,e ≃ T_int, the solutions are gained for a given value of T_int. Accordingly, these solutions should be viewed as an upper limit; disk truncation can be achieved when the accretion rate is lower than the dashed lines (the bottom left panel of Figure 2). Moreover, the negative and weak dependences on M_p arise due to a constant T_int for ${\dot{M}}_{{\rm{p}}}$ and L_acc, respectively. On the other hand, ${\dot{M}}_{{\rm{p}}}^{\mathrm{therm}}$ and ${L}_{\mathrm{acc}}^{\mathrm{therm}}$ are self-consistent solutions, and hence, the resulting values (denoted by the solid lines) are required to establish disk truncation due to magnetospheric accretion.

Thus, the planetary magnetic fields become strong enough to truncate circumplanetary disks when the disk accretion rate is lower than ${\dot{M}}_{{\rm{p}}}^{\mathrm{int}}$ when T_p,e ≃ T_int and when it becomes comparable to ${\dot{M}}_{{\rm{p}}}^{\mathrm{therm}}$ for T_p,e ≃ T_therm. It should be pointed out that these constraints are obtained as a function of M_p for the given values of R_in in this section. In the following sections, we use the solutions derived here to gain further constraints on R_in

Gas in magnetospheric flow is known to be heated to ∼10⁴ K for CTTSs, from which Hα is emitted (e.g., Muzerolle et al. 2001). While the origin of the heat source is still unclear (e.g., Hartmann et al. 2016), it probably ultimately stems from the magnetic fields of the accreting objects and/or accretion energy. For young giant planets, this heating should predominantly be attributed to the accretion energy when the disk-limited gas accretion results in high accretion rates; as discussed in Section 2.3, the energetics of the accretion processes can be self-regulated by the accompanying heating when T_p,e ≃ T_therm. We here consider this case and derive a constraint on R_in.

We first explore the properties of the gas in a magnetospheric flow. A magnetospheric flow carries the following flux of energy when the disks are truncated at R = R_in (e.g., Calvet & Gullbring 1998):

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{\dot{M}}_{{\rm{p}}}{v}_{\mathrm{sh}}^{2}=\displaystyle \frac{{{GM}}_{{\rm{p}}}{\dot{M}}_{{\rm{p}}}}{{R}_{{\rm{p}}}}\left(1-\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}\right),\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{v}_{\mathrm{sh}} & = & \sqrt{\displaystyle \frac{2{{GM}}_{{\rm{p}}}}{{R}_{{\rm{p}}}}\left(1-\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}\right)}\\ & \simeq & 1.3\times {10}^{2}{{\rm{km}}\,{\rm{s}}}^{-1}{\left(\displaystyle \frac{{f}_{{\rm{T}}}}{3/4}\right)}^{1/2}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/2}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-1/2}\end{array}\end{eqnarray} \tag{ 19 }$$

is the velocity of the accreted gas at the planetary surfaces, and

$$\begin{eqnarray}&&{f}_{{\rm{T}}}=1-\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}.\end{eqnarray} \tag{ 20 }$$

In the above equation, we adopt that R_p/R_in = 1/4, that is, f_T = 3/4. The sound speed of gas at planetary surfaces with T_p,e of a few 10³ K (see Figure 2) becomes a few km s⁻¹ and is much lower than v_sh. Therefore, shocks are produced at the planetary surfaces when the accreted gas arrives there. Using the strong-shock approximation, we write the shock temperature (T_sh) as

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{sh}} & = & \displaystyle \frac{3}{16}\displaystyle \frac{\mu {m}_{{\rm{H}}}}{{k}_{{\rm{B}}}}{v}_{\mathrm{sh}}^{2}\\ & \simeq & 4.0\times {10}^{5}{\rm{K}}\left(\displaystyle \frac{\mu }{1}\right)\left(\displaystyle \frac{{f}_{{\rm{T}}}}{3/4}\right)\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right),\end{array}\end{eqnarray} \tag{ 21 }$$

where μ is the mean molecular weight of the accreted gas, m_H is the mass of the hydrogen nucleons, and k_B is the Boltzmann constant. The value of μ varies from ∼0.53 to ∼1.28 for ionized to neutral gas at solar abundance. Previous studies confirm that T_sh is high enough to both dissociate molecular hydrogen and ionize atomic hydrogen, which leads to hydrogen line emission including Hα (e.g., Aoyama et al. 2020).

The generation of shocks at the planetary surfaces is very likely and may be the most plausible explanation for PDS 70 b/c, as discussed above. Nonetheless, it may be interesting to estimate the temperature of the gas in a magnetospheric flow and examine whether Hα emission is possible from the flow, as it is for CTTSs. One conservative estimate of the flow temperature may be obtained, assuming that the kinetic energy carried by the magnetospheric flow completely dissipates before shocks occur and that the emission would behave like a blackbody at all wavelengths. The resulting flow temperature (T_{flow, BB}) is computed as

$$\begin{eqnarray}&&4\pi {R}_{{\rm{p}}}^{2}{f}_{\mathrm{fill}}{\sigma }_{\mathrm{SB}}{T}_{\mathrm{flow},\ \mathrm{BB}}^{4}=\displaystyle \frac{1}{2}{\dot{M}}_{{\rm{p}}}{v}_{\mathrm{sh}}^{2},\end{eqnarray} \tag{ 22 }$$

where f_fill is the so-called filling factor and represents the fraction of the planet surface area at which the magnetospheric flow arrives from the inner edge of the circumplanetary disks. Since Equation (22) is a function of ${\dot{M}}_{{\rm{p}}}$, we consider both the cases that T_p,e ≃ T_int and T_p,e ≃ T_therm, as done in Section 2.4,

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{flow},\ \mathrm{BB}}^{\mathrm{int}} & \simeq & 3.5\times {10}^{3}{\rm{K}}{\left(\displaystyle \frac{{f}_{\mathrm{fill}}}{{10}^{-2}}\right)}^{-1/4}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-1/4}{\left(\displaystyle \frac{{f}_{{\rm{T}}}}{3/4}\right)}^{1/4}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{5/24}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-3/8}{\left(\displaystyle \frac{{T}_{\mathrm{int}}}{700{\rm{K}}}\right)}^{2/3}\\ & & \times {\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-7/8}\end{array}\end{eqnarray} \tag{ 23 }$$

when T_p,e ≃ T_int, where Equation (12) is used, and

$$\begin{eqnarray}&&\begin{array}{l}{T}_{\mathrm{flow},\ \mathrm{BB}}^{\mathrm{therm}}\simeq 6.4\times {10}^{3}{\rm{K}}{\left(\displaystyle \frac{{f}_{\mathrm{fill}}}{{10}^{-2}}\right)}^{-1/4}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-3/4}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{1/2}\\ \qquad \times \ {\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)}^{1/2}{\left(\displaystyle \frac{{f}_{{\rm{T}}}}{3/4}\right)}^{1/4}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{5/8}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-9/8}\\ \qquad \times \ {\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-21/8}\end{array}\end{eqnarray} \tag{ 24 }$$

when T_p,e ≃ T_therm, where Equation (13) is used.

One recognizes that T_{flow, BB} becomes lower than 10⁴ K for both cases. This implies that Hα emission from magnetospheric flow would be unlikely for accreting giant planets if the emission behaves like blackbody. The possibility of a deviation from a blackbody would be very likely, however; the column density of the accretion flow may not be high enough to achieve a blackbody radiation at all the wavelengths. Instead, the flow may only be optically thick for certain line emission, such as Hα. In fact, Hα emission is known to be optically thick at the gas number density of n_H > 10¹² cm⁻³ at a temperature of 8000 K for the accretion flow onto CTTSs. (e.g., Storey & Hummer 1995; Zhu 2015). We find that a similar situation would be possible for an accretion flow around accreting giants as

$$\begin{eqnarray}\begin{array}{rcl}{n}_{{\rm{H}}}^{\mathrm{int}} & = & \displaystyle \frac{{\dot{M}}_{{\rm{p}}}^{\mathrm{int}}}{4\pi {R}_{{\rm{p}}}^{2}{f}_{\mathrm{fill}}{m}_{{\rm{H}}}{v}_{\mathrm{sh}}}\\ & \simeq & 4.5\times {10}^{12}{\mathrm{cm}}^{-3}{\left(\displaystyle \frac{{f}_{\mathrm{fill}}}{{10}^{-2}}\right)}^{-1}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-1}{\left(\displaystyle \frac{{f}_{{\rm{T}}}}{3/4}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{-2/3}{\left(\displaystyle \frac{{T}_{\mathrm{int}}}{700{\rm{K}}}\right)}^{8/3}{\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-7/2}\end{array}\end{eqnarray} \tag{ 25 }$$

when T_p,e ≃ T_int, where Equation (12) is used, and

$$\begin{eqnarray}\begin{array}{rcl}{n}_{{\rm{H}}}^{\mathrm{therm}} & \simeq & 4.7\times {10}^{13}{\mathrm{cm}}^{-3}{\left(\displaystyle \frac{{f}_{\mathrm{fill}}}{{10}^{-2}}\right)}^{-1}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-3}\\ & & \times {\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{2}{\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)}^{2}{\left(\displaystyle \frac{{f}_{{\rm{T}}}}{3/4}\right)}^{-1/2}\\ & & \times \left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right){\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-3}{\left(\displaystyle \frac{{R}_{\mathrm{in}}/{R}_{{\rm{p}}}}{4}\right)}^{-21/2}\end{array}\end{eqnarray} \tag{ 26 }$$

when T_p,e ≃ T_therm, where Equation (13) is used.

We now turn our attention to deriving a constraint on R_in. As discussed in Section 2.3, disk-limited gas accretion becomes energetic enough to self-regulate the thermal properties of the other processes taking place in this stage when T_p,e ≃ T_therm. For this case, the conservation of energy dictates

$$\begin{eqnarray}&&{L}_{\mathrm{gas}}\geqslant \displaystyle \frac{1}{2}{\dot{M}}_{{\rm{p}}}{v}_{\mathrm{sh}}^{2},\end{eqnarray} \tag{ 27 }$$

where the inequality sign appears as some energy may radiate away from magnetospheric flow (see Section 2.1). The above equation is rewritten as (see Equations (2) and (18))

$$\begin{eqnarray}&&1-\displaystyle \frac{3}{2}\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}\right)\left[1-\displaystyle \frac{2}{3}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}\right)}^{1/2}\right]\geqslant \left(1-\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{\mathrm{in}}}\right).\end{eqnarray} \tag{ 28 }$$

This reads R_in ≤ 4R_p. Note that the above condition is only applicable to the case that T_p,e ≃ T_therm; when T_p,e ≃ T_int, the planetary magnetic fields, and hence, the disk truncation radius, are determined by the intrinsic temperature (Equations (8) and (31), also see Section 2.6). These quantities are controlled by earlier formation histories and may provide additional heating for the accretion flow.

In summary, a magnetospheric flow around accreting giant planets can result in Hα emission either via an accretion shock at the planetary surfaces or via an accretion flow that would be heated by inefficient cooling of certain lines. Due to the conservation of energy, R_in ≤ 4R_p when disk-limited gas accretion self-regulates the energetics of the processes operating at this stage.

We have focused on gas accretion flow in the vicinity of magnetized giant planets so far. We here consider a more global configuration of the accretion flow. In particular, we explore what the gas flow from parental circumstellar disks onto accreting planets and the surrounding circumplanetary disks looks like. This consideration becomes important when planets and their circumplanetary disks are embedded in the circumstellar disks. For this case, we can obtain another constraint on R_in that comes from the surrounding environment (e.g., stellar accretion rates).

The gas accretion flow from circumstellar disks onto circumplanetary disks and/or the host planets is currently poorly constrained. The primary reason is that observations of circumplanetary disks are so far very limited; the disks around PDS 70 b/c are the only examples so far (e.g., Isella et al. 2019; Benisty et al. 2021). The lack of observations hinders a specification of the disk properties, and hence, a development of reliable models. Under this circumstance, we only consider the overall structure of the gas accretion flow in this work.

We make use of the approach of Tanigawa & Tanaka (2016) to estimate how much of the gas is delivered from the parental circumstellar disks to the system of planets and their circumplanetary disks. In the approach, the results of two different hydrodynamical simulations are coupled together; one type of simulations derives a formula for the accretion rate onto a system like this (Tanigawa & Watanabe 2002), and the other simulations compute the reduction factor of the surface density of circumstellar disks, which is caused by disk–planet interaction (Kanagawa et al. 2015). The resulting gas flow rate (${\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}$) is given as

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}=\displaystyle \frac{8.5}{3\pi }\left(\displaystyle \frac{{c}_{{\rm{s}}}^{\mathrm{CSD}}}{{v}_{\mathrm{Kep}}^{\mathrm{CSD}}}\right){\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{s}}}}\right)}^{-2/3}{\dot{M}}_{{\rm{s}}},\end{eqnarray} \tag{ 29 }$$

where ${c}_{{\rm{s}}}^{\mathrm{CSD}}$ and ${v}_{\mathrm{Kep}}^{\mathrm{CSD}}$ are the sound speed and the Keplerian velocity of the circumstellar disk gas at the position of planets, and M_s and ${\dot{M}}_{{\rm{s}}}$ are the mass of the central star and the disk accretion rate onto the star, respectively.

The value of ${\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}$ becomes comparable to the accretion rate onto planets (${\dot{M}}_{{\rm{p}}}$) when the accretion flow is in a steady state. There is no guarantee that accreting giant planets achieve this state. However, Hasegawa et al. (2021) have recently found that this might be the case for PDS 70 b/c; their calculations show that

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}} & \simeq & 1.8\times {10}^{-7}{M}_{{\rm{J}}}\,{{\rm{yr}}}^{-1}\left(\displaystyle \frac{{c}_{{\rm{s}}}^{\mathrm{CSD}}/{v}_{\mathrm{Kep}}^{\mathrm{CSD}}}{8.9\times {10}^{-2}}\right)\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{-2/3}\left(\displaystyle \frac{{\dot{M}}_{{\rm{s}}}}{1.4\times {10}^{-10}{M}_{\odot }\,{{\rm{yr}}}^{-1}}\right),\end{array}\end{eqnarray} \tag{ 30 }$$

where ${c}_{{\rm{s}}}^{\mathrm{CSD}}/{v}_{\mathrm{Kep}}^{\mathrm{CSD}}=8.9\times {10}^{-2}$ and M_s = 0.76 M_⊙ are adopted, following Keppler et al. (2018), who simulate the properties of the circumstellar disk around PDS 70, and the value of ${\dot{M}}_{{\rm{s}}}$ is taken from Thanathibodee et al. (2020), suggesting that ${\dot{M}}_{{\rm{s}}}$ of PDS 70 lies within the range of 0.6–2.2 × 10⁻¹⁰ M_⊙ yr⁻¹.

It is noticeable that within the range of ${\dot{M}}_{{\rm{s}}}$, the resulting value of ${\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}$ is comparable to the value of ${\dot{M}}_{{\rm{p}}}\simeq {10}^{-8}\mbox{--}{10}^{-7}{M}_{{\rm{J}}}\,{{\rm{yr}}}^{-1}$, which is estimated from the observed Hα emission for PDS 70 b/c (Hasegawa et al. 2021 and references herein). This implies that the steady-state accretion assumption may not be unreasonable for accreting giant planets, at least during certain formation stages.

Motivated by the finding, we use the assumption and derive a constraint on R_in. To proceed, we consider two limiting cases: T_p,e ≃ T_int and T_p,e ≃ T_therm, and we equate ${\dot{M}}_{{\rm{p}}}^{\mathrm{int}}$ and ${\dot{M}}_{{\rm{p}}}^{\mathrm{therm}}$ with ${\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}$. When T_p,e ≃ T_int,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{R}_{\mathrm{in}}^{\mathrm{int}}}{{R}_{{\rm{p}}}} & \simeq & 2.0{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-2/7}{\left(\displaystyle \frac{{h}_{0}}{0.05}\right)}^{-2/7}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/7}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{3/7}{\left(\displaystyle \frac{{T}_{\mathrm{int}}}{700{\rm{K}}}\right)}^{16/21}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{s}}}}{1{M}_{\odot }}\right)}^{-4/21}{\left(\displaystyle \frac{{\dot{M}}_{{\rm{s}}}}{{10}^{-10}{M}_{\odot }\,{{\rm{yr}}}^{-1}}\right)}^{-2/7}\\ & & \times {\left(\displaystyle \frac{{r}_{{\rm{p}}}}{50\,{\rm{au}}}\right)}^{-1/14},\end{array}\end{eqnarray} \tag{ 31 }$$

where Equations (12) and (29) are used, and when T_p,e ≃ T_therm,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{R}_{\mathrm{in}}^{\mathrm{therm}}}{{R}_{{\rm{p}}}} & \simeq & 4.0{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-2/7}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{4/21}{\left(\displaystyle \frac{{h}_{0}}{0.05}\right)}^{-2/21}\\ & & \times {\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)}^{4/21}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{13/63}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-1/7}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{s}}}}{1{M}_{\odot }}\right)}^{-4/63}{\left(\displaystyle \frac{{\dot{M}}_{{\rm{s}}}}{{10}^{-10}{M}_{\odot }\,{{\rm{yr}}}^{-1}}\right)}^{-2/21}\\ & & \times {\left(\displaystyle \frac{{r}_{{\rm{p}}}}{50\,{\rm{au}}}\right)}^{-1/42},\end{array}\end{eqnarray} \tag{ 32 }$$

where Equations (13) and (29) are used. Note that in the above equations, ${c}_{{\rm{s}}}^{\mathrm{CSD}}/{v}_{\mathrm{Kep}}^{\mathrm{CSD}}\equiv {h}_{0}{({r}_{{\rm{p}}}/1\,{\rm{au}})}^{1/4}$ is assumed, where h₀ = 0.05, and r_p is the position of planets. Moreover, the dependence on 1 − f_in is weak, and an intermediate value of 3/4 is used in Equation (32); the value of 1 − f_in varies from 1/2 to 1.

Figure 3 visualizes the change in ${R}_{\mathrm{in}}^{\mathrm{int}}$ and ${R}_{\mathrm{in}}^{\mathrm{therm}}$ as a function of M_p for a given value of ${\dot{M}}_{{\rm{s}}};$ since the dependence on other parameters including r_p is very weak (see Equations (31) and (32)), we focus on M_p and ${\dot{M}}_{{\rm{s}}}$. Furthermore, 1 − f_in is set at 3/4, as was done in Equation (32). The resulting trends can be understood readily; when the stellar accretion rates are high, the ram pressure becomes strong, and hence, the inner disk edge is located close to the host planets. For the M_p dependence, a monotonic increase of R_in arises due to disk–planet interaction (Equation (29)); massive planets open up a deep gap in their parental circumstellar disks, which decreases the gas flow (${\dot{M}}_{{\rm{p}}}^{\mathrm{CDS}}$) onto these planets and the surrounding circumplanetary disks. As a result, R_in expands due to the low ram pressure. Our calculations show that when T_p,e ≃ T_therm, $1\lt {R}_{\mathrm{in}}^{\mathrm{therm}}\leqslant 4$ for planets with M_p ≤ 10M_J, suggesting that magnetospheric accretion is possible for a wide range of parameters. On the other hand, when T_p,e ≃ T_int, ${R}_{\mathrm{in}}^{\mathrm{int}}\lt 1$ for high stellar accretion rates. Therefore, magnetospheric accretion only occurs in the later stage of disk evolution.

Zoom In Zoom Out Reset image size

The condition (${R}_{\mathrm{in}}^{\mathrm{int}}/{R}_{{\rm{p}}}\gt 1$) needed for the case that T_p,e ≃ T_int leads to a constraint on ${\dot{M}}_{{\rm{s}}}$ as

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{{\rm{s}}}\lt 1.1 & & \times {10}^{-9}{M}_{\odot }\,{{\rm{yr}}}^{-1}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-1}{\left(\displaystyle \frac{{h}_{0}}{0.05}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/2}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{3/2}{\left(\displaystyle \frac{{T}_{\mathrm{int}}}{700{\rm{K}}}\right)}^{8/3}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{s}}}}{1{M}_{\odot }}\right)}^{-2/3}{\left(\displaystyle \frac{{r}_{{\rm{p}}}}{50\,{\rm{au}}}\right)}^{-1/4},\end{array}\end{eqnarray} \tag{ 33 }$$

where Equations (12) and (29) are used. This is equivalent to the condition that ${\dot{M}}_{{\rm{p}}}^{\mathrm{int}}\gt {\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}$ at R_in = R_p. When T_p,e ≃ T_therm, the required condition (${R}_{\mathrm{in}}^{\mathrm{therm}}/{R}_{{\rm{p}}}\leqslant 4$) is rewritten as

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{{\rm{s}}} & \geqslant & 9.0\times {10}^{-11}{M}_{\odot }\,{{\rm{yr}}}^{-1}{\left(\displaystyle \frac{{f}_{\mathrm{ram}}}{1/\sqrt{2}}\right)}^{-3}{\left(\displaystyle \frac{{f}_{{\rm{L}}}}{1/2}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{{h}_{0}}{0.05}\right)}^{-1}{\left(\displaystyle \frac{1-{f}_{\mathrm{in}}}{3/4}\right)}^{2}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{13/6}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-3/2}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{s}}}}{1{M}_{\odot }}\right)}^{-2/3}{\left(\displaystyle \frac{{r}_{{\rm{p}}}}{50\,{\rm{au}}}\right)}^{-1/4},\end{array}\end{eqnarray} \tag{ 34 }$$

where Equations (13) and (29) are used. This can also be obtained from the condition that ${\dot{M}}_{{\rm{p}}}^{\mathrm{therm}}\leqslant {\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}$ at R_in = 4R_p.

In the following section, we use Equations (33) and (34) and predict when hydrogen lines can be emitted from young giant planets via magnetospheric accretion, either due to an accretion shock or to the inefficiently cooled accretion flow.

Armed with the equations derived in the above sections, we are now ready to explore the line luminosity of the hydrogen emission originating from giant planets undergoing magnetospheric accretion. In order to compute the line luminosity, we strongly rely on the relationships between the line and accretion luminosities that are obtained by previous studies: Aoyama et al. (2021) for emission from an accretion shock, and Alcalá et al. (2017) for emission from an accretion flow along magnetospheres. The former computed the relationship theoretically, and the latter obtained it observationally from CTTSs.

We first summarize the key equations for computing the accretion luminosity. Since we target planets surrounded by their circumplanetary disks, which are embedded in their parental circumstellar disks, we assume that $\dot{{M}_{{\rm{p}}}}\simeq {\dot{M}}_{{\rm{p}}}^{\mathrm{CSD}}$, as discussed in Section 2.6. Then, the accretion luminosity is written as

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{acc}} & \simeq & 5.2\times {10}^{-4}{L}_{\odot }\left(\displaystyle \frac{{h}_{0}}{0.05}\right){\left(\displaystyle \frac{{M}_{{\rm{p}}}}{10{M}_{{\rm{J}}}}\right)}^{1/3}{\left(\displaystyle \frac{{R}_{{\rm{p}}}}{1.5{R}_{{\rm{J}}}}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{s}}}}{1{M}_{\odot }}\right)}^{2/3}\left(\displaystyle \frac{{\dot{M}}_{{\rm{s}}}}{{10}^{-10}{M}_{\odot }\,{{\rm{yr}}}^{-1}}\right){\left(\displaystyle \frac{{r}_{{\rm{p}}}}{50\,{\rm{au}}}\right)}^{1/4},\end{array}\end{eqnarray} \tag{ 35 }$$

where Equations (1) and (29) are used.

We then consider two limiting cases for T_p,e: when planets undergo magnetospheric accretion and their effective temperature is given as T_p,e ≃ T_int, the accretion luminosity can only be emitted when $1\lt {R}_{\mathrm{in}}^{\mathrm{int}}/{R}_{{\rm{p}}}$ (see Equations (31) and (33)). The magnetic fields of these planets are written by Equation (8). On the other hand, when planets with an effective temperature of T_p,e ≃ T_therm (see Equation (15)) experience magnetospheric accretion, then the accretion luminosity can only be emitted when $(1\lt ){R}_{\mathrm{in}}^{\mathrm{therm}}/{R}_{{\rm{p}}}\leqslant 4$ (see Equations (32) and (34)). Their magnetic fields are given by Equation (16).

Separating the contributions of T_p,e (T_int versus T_therm in Equation (4)) allows us to identify a parameter space wherein each contribution becomes dominant. This is clearly shown in Figure 4, which plots the conditions of which the value of L_acc can be emitted from giant planets via magnetospheric accretion. The value of L_acc increases monotonically with increasing M_p and ${\dot{M}}_{{\rm{s}}}$, which is obvious from Equation (35). As anticipated from the above discussion, the ${M}_{{\rm{p}}}-{\dot{M}}_{{\rm{s}}}$ parameter space is divided into three regions: the region above the dashed line, where disk-limited gas accretion is energetic enough to self-regulate the resulting line emission (i.e., T_p,e ≃ T_therm), the region below the solid line, where early formation histories play an important role for L_acc even at the disk-limited gas accretion stage (i.e., T_p,e ≃ T_int), and the intermediate region, where both the cases are possible. Our calculations show that magnetospheric accretion leads to L_acc that can be high enough to be observed for certain combinations of parameters (i.e., high M_p and high ${\dot{M}}_{{\rm{s}}}$).

Zoom In Zoom Out Reset image size

We now compute the line luminosity (L_line) of hydrogen emission using the following equation:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}({L}_{\mathrm{acc}}/{L}_{\odot })=a\times {\mathrm{log}}_{10}({L}_{\mathrm{line}}/{L}_{\odot })+b,\end{eqnarray} \tag{ 36 }$$

where the fitting parameters (a and b) are summarized in Table 2. As an example, we consider three lines that tend to be observed readily (see Table 2). Similar calculations are straightforward for other lines (see table 1 of Aoyama et al. 2021, where the values of a and b are tabulated for other lines).

Figure 5 shows the resulting line luminosities in the ${M}_{{\rm{p}}}-{\dot{M}}_{{\rm{s}}}$ parameter space. The value of r_p = 50 au is chosen, as in Figure 4. It is obvious from Table 2 that the Hα line luminosity is dimmer by more than one order of magnitude than the accretion luminosity, and the line luminosity is higher for the accretion flow case than the accretion shock case. For Paβ and Brγ, similar trends are confirmed, while L_Paβ and L_Brγ are lower by more than two and three orders of magnitude than L_acc, respectively.

Zoom In Zoom Out Reset image size

In summary, one can predict the line luminosities of hydrogen that will be emitted from young giant planets due to magnetospheric accretion when the accretion rate onto the host star and the mass and position of the planets are estimated. The observed value of the line luminosities (and the line ratios) can be used as a diagnostics to identify where the emission originates (planetary surface versus accretion flow) and how the emission is produced (accretion shock versus accretion heating). The specification of the stellar accretion rates and planet properties enables a determination of the ultimate cause of why these planets undergo magnetospheric accretion, namely, the magnetism of young accreting giant planets; disk-limited gas accretion is energetic enough to trigger it, or early planet formation processes keep planets hot enough.

We observed HD 163296 on 2021 May 8 UT with Subaru/SCExAO+VAMPIRES under the NASA-Keck time exchange program (PID 61/2021A_N200: PI—Hasegawa).

VAMPIRES has two detectors that can take different images with two filters simultaneously, and it is capable of mitigating aberrations between the detectors by switching the filters (double-differential calibration; Norris et al. 2015). When conducting Hα imaging with VAMPIRES, we used narrow-band filters for Hα (λ_c = 656.3 nm, Δλ = 1.0 nm) and the adjacent continuum (λ_c = 647.68 nm, Δλ = 2.05 nm), which allowed us to effectively subtract the continuum components from the Hα image (spectral differential imaging, SDI; Smith 1987). During the observations, we repeated two states for the double-differential calibration; in State 1, cam1 is used for the continuum and cam2 for Hα, and in State 2, the setup is reversed. Note that due to an instrumental constraint, we used a smaller field of view (FoV: ∼15 × 15) than the maximum FoV of VAMPIRES (∼3'' × 3'').

The single-exposure time was 50 ms in the first sequence (9 cubes in both State 1 and State 2), and from the second sequence, we changed it to 40 ms in order to avoid saturation. The difference of the single-exposure time was corrected before postprocessing. The total integration time corresponds to 2806 s and 2369 s, and the field rotation angle gains are ∼64° and 62° for States 1 and 2, respectively.

The VAMPIRES data format is a cube consisting of an image and short exposures (2001 exposures per cube). We first subtracted dark from each exposure and conducted a point-spread function (PSF) fitting of the continuum frames by a 2D Gaussian for the frame selection. The typical full width at half maximum (FWHM) of the PSF was measured at ∼20–25 mas with a pixel scale of 6.24 ± 0.01 mas pix⁻¹ (Currie et al. 2022). We investigated the fitted peaks and removed a few data cubes that did not exhibit the typical peaks due to poor adaptive optics corrections. We then empirically selected the 80 percentile of the fitted peak values in a cube and then combined the selected exposures into an image after aligning the centroid of the PSFs (see Figure 2 of Uyama et al. 2020).

In order to remove the stellar halo and to search for faint accretion signatures by postprocessing, we followed the postprocessing methods of Uyama et al. (2022; see their Section 2.1 for details and references herein), using angular differential imaging (ADI; Marois et al. 2008), SDI with the two filters, and the VAMPIRES double-differential calibration techniques. We scaled the continuum images by calculating a scaling factor from the comparison between the photometry (aperture radius=10 FWHM) of the Hα and continuum filters and by correcting the wavelengths so that they were appropriate for the reference PSF of the SDI reduction. In the ADI reduction, we used pyklip packages (Wang et al. 2015), which make a reference PSF by Karhunen-Loève image projection (KLIP; Soummer et al. 2012). Note that we adopted an aggressive ADI reduction to explore the faintest possible accretion signatures, and this setting can attenuate extended features. Therefore, we do not discuss the Hα jet features that are present in Xie et al. (2021). We also note that the A4 knot detected by Xie et al. (2021) is beyond the FoV in our observations. We then applied SDI using the ADI residuals of Hα and the scaled continuum at States 1 and 2, and we finally performed a double-differential calibration.

We did not find any companion candidates within 07. ¹⁵ This is clearly shown in Figure 6. We calculated standard deviations within annular regions after convolving the output image with a radius of FWHM/2. The image is then compared with photometry of the central star with an aperture radius with a FWHM/2 at the Hα filter for a contrast limit (Figure 7). We also took into account throughput loss caused by the ADI reduction, where false PSFs are injected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to directly compare the observation results with theoretical predictions, we convert the contrast limit into the Hα flux limit. The conversion is done by referring to the continuum flux of the central star and taking the Hα/continuum ratio from our observations into account. The resulting conversion factor is ∼1.3 × 10⁻¹⁰. Note that the detection limit corresponds to the integrated line flux of Hα as our observations cannot resolve the line, and thus we do not take the line profile into account for the comparison between the observational results and our model. The VAMPIRES Hα narrow-band filter with 1.0 nm corresponds to a velocity coverage of±100 km s⁻¹, which is well above the possible maximum gas velocity around a Jovian protoplanet (see also Appendix A in Uyama et al. 2020). Sitko et al. (2008) reported that HD 163296 exhibits no variability within 10%, except for a significant variability in the near-IR wavelengths per 16 yr. We therefore used the Gaia G-band flux (4.81 erg s⁻¹ cm⁻² μm⁻¹; Gaia Collaboration et al. 2018) as the continuum. The Hα/continuum ratio is estimated at 2.66 by comparing the photometry of the central star between the Hα and the continuum filters, and assuming the multiplied value as the HD 163296 flux at Hα. An aperture radius of 10 FWHM is used in the above conversion.

Figure 7 shows the 5σ detection limit as a function of the distance from the host star. Note that the detection limit is not computed by Xie et al. (2020, 2021), where the MUSE data taken toward HD 163296 are analyzed; these data include instrumental noises, which makes it hard to accurately estimate the detection limit. No direct comparison between our detection limit and the MUSE limit is therefore made in this work.

In the following sections, we use the above detection limit and apply our theoretical predictions to the HD 163296 system.

Before comparing our predictions with the observational results, we here consider the effect of extinction.

Extinction occurs when gas and/or dust are present between the emitting sources and the observer, which can potentially reduce the observed line flux significantly from the intrinsically emitted flux. Its value (A_λ ) measured in magnitude at a wavelength λ is defined as (e.g., Draine 2011)

$$\begin{eqnarray}&&{A}_{\lambda }\equiv 2.5{\mathrm{log}}_{10}\left({F}_{\lambda }/{F}_{\lambda }^{\mathrm{obs}}\right),\end{eqnarray} \tag{ 37 }$$

where ${F}_{\lambda }^{\mathrm{obs}}$ is the actually observed flux, and F_λ is the intrinsic flux emitted from the sources before extinction comes into play. In this work, F_λ corresponds to the theoretically computed value.

The value of A_λ is quantified relatively well for star-forming environments (e.g., Draine 2011). In fact, A_λ is written as

$$\begin{eqnarray}&&{A}_{\lambda }={N}_{{\rm{H}}}/{K}_{\lambda },\end{eqnarray} \tag{ 38 }$$

where N_H is the total column density of hydrogen distributed between the sources and the observer, and K_λ is the conversion coefficient. For a diffuse interstellar medium (ISM) and molecular clouds, the value of K_V is known to be on the order of 10²¹ mag⁻¹ cm⁻² at a visual wavelength (i.e., λ = 0.55 μm; e.g., Bohlin et al. 1978; Olofsson & Olofsson 2010). On the other hand, the effect of extinction is poorly constrained for young accreting giant planets; observations of these planets are currently very rare, and hence, the emitting environment remains to be studied. Theoretically, extinction originating from gas is expected to be small at least at Hα (e.g., Marleau et al. 2022). However, the dust opacity can be non-negligible at optical and infrared wavelengths (e.g., Sanchis et al. 2020). In this work, we therefore attempt to compute the value of K_λ using the Hα observations obtained for PDS 70 b/c and our theoretical models.

Table 3 summarizes the input parameters and computed quantities. We use Equations (35) and (36) to compute (the theoretically predicted) intrinsic line flux. The values of extinction and the coefficient are then calculated from Equations (37) and (38), respectively. The (observed) input parameters are taken from Hashimoto et al. (2020). In addition, the stellar mass, the mass of planets b and c, and the surface density of the circumstellar disk around the planet positions are assumed to be M_s = 0.85M_⊙, M_p ∼ 2M_J, and ${{\rm{\Sigma }}}_{{\rm{d}}}^{\mathrm{CSD}}\,\sim 0.1$ g cm⁻², respectively, following Keppler et al. (2019). The last quantity is used to compute ${N}_{{\rm{H}}}(={{\rm{\Sigma }}}_{{\rm{d}}}^{\mathrm{CSD}}/{m}_{{\rm{H}}})$. The distance of PDS 70 from Earth is set at 113 pc (Hashimoto et al. 2020).

Our calculations show that even when the wavelength dependence of A_λ ( ∝ λ^−1.75) is taken into account (e.g., Draine 2011), the coefficient K_Hα is about a few times higher than the value obtained in star-forming environments. This is likely to be reasonable as gas contributing to extinction for accreting giants may come from the surface layer of the parental circumstellar disks; this disk gas may be poor in the dust abundance due to dust settling and growth compared with the ISM gas. The value of extinction itself can nonetheless be higher than that of star-forming environments simply because N_H may be much higher in planet-forming environments. It should be pointed out that our estimate of A_λ for the accretion shock case is comparable to that of Hashimoto et al. (2020), where the extinction values are derived from the line flux ratio of the observed Hα and nondetected Hβ.

In the following section, we use the computed value of K_Hα to take into account extinction for the HD 163296 system.

We finally compare our theoretical predictions made in Section 2 with our observational results obtained toward the HD 163296 system.

Figure 8 shows the results. The line flux of Hα for the planet candidates is computed using Equations (35), (36), (37), and (38). The properties of these candidates are summarized in Table 1. The value of ${{\rm{\Sigma }}}_{{\rm{d}}}^{\mathrm{CSD}}\sim 0.5$ g cm⁻² is used, following Isella et al. (2016). The error bars come from the ranges of planet mass and positions and from the variation in K_Hα (see Table 3). Our calculations show that the observational results are not sensitive enough to reliably examine the theoretical predictions developed in Section 2; an examination like this requires that the observational sensitivity should be increased by one order of magnitude or more. We also find that an accretion shock leads to a higher observed flux than an accretion flow, which is expected from the value of K_Hα (Table 3). In addition, we have confirmed that the emission resides in the region where T_p,e ≃ T_therm (see Figure 5), and hence, if Hα were observed toward the HD 163296 system, then the line could be used as a direct probe of the disk-limited gas accretion stage of giant planet formation.

Zoom In Zoom Out Reset image size

It can thus be concluded that the hydrogen emission lines, especially Hα, are useful tracers of whether giant planets undergo magnetospheric accretion in their final formation stages. However, the current observational capability may not be high enough to reliably test the emission mechanisms (e.g., accretion shock versus accretion flow); when the planets are embedded in actively accreting circumstellar disks, the emission from the planets itself can be strong. Since the emission is the direct outcome of a high accretion flow onto the planets, the flow in turn attenuates the observed flux significantly. When planets are located in disks with low stellar accretion rates, the emission becomes weaker, which simply makes it difficult to observe. An improvement of the observational sensitivity by a factor of 10 or more will open up a promising window to carefully investigate the final giant planet formation stage.

As described above, the Hα lines tend to be significantly extincted. We here explore other lines (e.g., Paβ and Brγ) that are less strongly attenuated by the magnetospheric accretion flow.

Figure 9 shows the resulting line flux for Paβ and Brγ. In these calculations, we adopt the same input parameters as in Section 3.4. We find that the observable line flux for Paβ and Brγ should be much higher than that of Hα. This arises because extinction is a decreasing function of λ (i.e., A_λ ∝ λ^−1.75) in our calculations. We have confirmed that contamination from the continuum emission by the planets and the disk is negligible for the HD 163296 system.

Zoom In Zoom Out Reset image size

Thus, the intrinsic hydrogen emission lines originating from accreting giants are weaker with increasing wavelengths (see Figure 5). However, the effect of extinction also becomes weaker for longer wavelengths. As a result, the observability of these lines (e.g., Paβ and Brγ) becomes higher than that of Hα. Multiband observations will expand the possibility of discovering and characterizing young giant planets that are embedded in the parental circumstellar disks that may undergo magnetospheric accretion.